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Numerical simulations of assemblies of grains under cyclic loading exhibit "granular ratcheting" : 
a small net deformation occurs with each cycle, leading to a linear accumulation of deformation with 
cycle number. We show that this is due to a curious property of the most frequently used models 
of the particle-particle interaction: namely, that the potential energy stored in contacts is path- 
dependent. There exist closed paths that change the stored energy, even if the particles remain 
in contact and do not slide. An alternative method for calculating the tangential force removes 
grarmlar ratcheting. 

PACS numbers: 45.70.-n,45.10.-b,83.10.Rs,45.20.dh 
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I. INTRODUCTION 

Granular ratcheting refers to the slow, linear accumu- 
lation of strain in a granular sample under cyclic loading. 
Several versions of this phenomena have been identified. 
The first variant to be found occurs when the loaded sam- 
ple reaches the critical state once per cycle. The mech- 
anism is easily understood: the material flows while it 
is in the critical state, giving rise to a deformation that 
accumulates with cycle number. However, ratcheting can 
also appear even when the sample never reaches the crit- 
ical state [H . In the following, we discuss exclusively this 
second type of ratcheting. 

Ratcheting in the absence of a critical state has also 
been observed in numerical simulations [3, 0j S ■ This 
is a very promising development, for one has access to all 
the quantities in numerical simulations, and it is usually 
possible to identify the origin of the phenomena. Once 
this has been done, one can then ask if the cause of the 
phenomena in the simulations is related to the cause in 
the experiments. 

Numerical studies have already provided many clues to 
granular ratcheting. The important role of sliding con- 
tacts has been pointed out granular ratcheting has 
been delimited from other possible behaviors Q , and the 
influence of various parameters has been studied 0, Q . 
One finding of these studies is that granular ratcheting 
is a quasi-static phenomena. Specifically, if one lets the 
frequency of the oscillating force tend to zero while keep- 
ing all other parameters the same, the deformation per 
cycle approaches a constant. 

What is still missing is an understanding of granular 
ratcheting on the micro-mechanical level: How exactly 
does the phenomenon arise from interaction of individ- 
ual particles? Is it possible to modify the particle inter- 
action law to eliminate ratcheting? What is the simplest 
system needed to produce ratcheting? We answer these 
questions in this paper. 

In Sec. in] we show that ratcheting can be obtained 
with only 16 particles. Ratcheting occurs if a single con- 



tact becomes sliding. Previously, ratcheting was linked 
to the presence of sliding contacts, but this is the first 
time that it is shown that only a single contact is nec- 
essary. We also show that the application of unconven- 
tional boundary conditions can lead to ratcheting even 
when there are no sliding contacts. In Sec. IIIH we show 
how granular ratcheting arises from the way tangential 
forces are calculated. In Sec. lIV Al we present an alterna- 
tive method that does not exhibit ratcheting. Finally, in 
the appendix, we show how stiffness matrix theory can 
illuminate some aspects of the problem. 



II. DESCRIPTION OF GRANULAR 
RATCHETING 

A. Model definition 

In this section, we present a very brief description of 
granular ratcheting, since more complete discussions al- 
ready exist 0] . Granular ratcheting is observed in biaxial 
or triaxial tests, where a granular sample is enclosed in 
a test chamber, and subjected to a uniform pressure and 
a cyclic load. We consider here exclusively the two di- 
mensional version of these experiments, often called the 
"biaxial box" , where a granular sample composed of disks 
is enclosed in a rectangular box of dimensions x Ly, 
with forces and Fy exerted on the walls. The forces 
are 



F^ 



PoLy, 



Fy = [Po + qit)] 



(1) 



where Pq is the pressure exerted on the sample, and q(t) is 
a periodic function, usually sinusoidal. In the simulations 
presented here, q{t) — Acr(l — cos ujt). One usually uses 
deviatoric strain 
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Lx 



(2) 



to characterize the deformation. (Here Lxo and Lyo are 
the lengths of the system at the beginning of the Simula- 
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tion.) 

We use a common numerical model of granular materi- 
als: grains are represented by disks who repel each other 
when they overlap. Thus whenever two disks touch each 
other, they exert a repulsive force Fn at the point of con- 
tact, directed normal to the particle surfaces. Fn is an 
increasing function of the overlap If the surfaces 

of two touching disks move relative to each other in the 
tangential direction, a second force Ft arises, directed 
tangent to the particle surfaces. F„ and Ft are called the 
normal and tangential components of the contact force. 
In addition to these forces, some damping forces are in- 
cluded to remove energy injected by the loading. 

The contact force is subjected to two constraints, 
namely 



Fn>0, ^iFn>\Ft\ 



(3) 



The first condition excludes cohesion, and the second is 
the Coulomb friction law. The constant /i is the Coulomb 
friction coefhcient. Contacts where |Ft| = jiFn are called 
sliding contacts, and those where \Ft\ < /iF„ are called 
non-sliding. 

All studies of granular ratcheting use this model, ex- 
cept sometimes polygons are used instead of disks [3j . 



B. Ratcheting with sixteen particles 

If one wishes to approximate the continuum-like behav- 
ior of soils, simulations with large numbers of particles 
are necessary. Therefore, granular ratcheting has been 
studied in assemblies of hundreds or thousands of parti- 
cles. In this paper, however, we wish simply to discover 
the origin of the phenomena, so it is useful to consider 
small numbers of particles. In this section, we study an 
assembly of sixteen particles that exhibits granular ratch- 
eting. The normal force is taken to be proportional to 
the overlap area, as in Ref. [5|. 

In Fig. [T] we show a plot of 7 vs q for a biaxial test 
performed on sixteen circular particles. In the first cycle, 
7 increases to about 0.008. During subsequent cycles, 
the system appears to trace out a four-sided polygon in 
the q-j plane. However, the path is not quite a polygon, 
because the system does not quite return to its starting 
point after one cycle, but to one where 7 is slightly larger. 
This is made obvious in Fig. [21 where 7 is plotted at 
t — nT, were T is the period of the cyclic loading, and 
n = 0,1,2 .. .. A small, linear increase of 7 with cycle 
number is visible. This is granular ratcheting. 

During one cycle, all the contacts remain non-sliding, 
except for one, which becomes sliding twice per cycle. 
This single contact is responsible for granular ratcheting, 
for if we inhibit sliding at this contact by increasing fi, 
granular ratcheting stops. However, we show below that 
ratcheting can also occur without sliding contacts with 
slightly different boundary conditions. 

In Fig. [HI we show the trajectory of the sliding contact 
in its {Fn,Ft) plane. The equalities \Ft\ = /xF„ are also 
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FIG. 1: Stress-strain curve of the system with 16 particles 
discussed in Sec. Ill Bl 
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FIG. 2: Accumulated strain for the system shown in Fig. [1] 



shown on the graph, and form a cone, with the vertex at 
the origin. The conditions in Eq. ^ mean that {Fn,Ft) 
must always lie within this cone. As one can see, the 
ratcheting contact's trajectory is a trapezoid, with the 
four corners labeled A, B, C , and D. The two paral- 
lel line segments correspond to the change in force when 
all contacts are non-sliding. Line segments BC and DA 
lie on the sides of the cone \Ft\ = ^Fn, and correspond 
to times when the contact is sliding with Ft = fJ.F„ or 
Ft = —iJ-Fn respectively. The trajectory is not quite a 
trapezoid, because after one cycle, the point does not re- 
turn to A, but arrives at A' , a bit closer to the origin 
than, but very close to, A. After the following cycle, the 
system has again shifted towards the origin by the same 
amount. This shift has its origin in the tiny displace- 
ment that occurs with each cycle - the sliding contact is 
gradually opening. 

The non-sliding contacts in the packing also trace out 
trapezoids, but their edges do not intersect the cone 
\Ft\ = nFn- Some examples are shown in Fig. [H The 
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FIG. 3: A sketch of the shding contact's trajectory in its 
{Fn,Ft) plane. The diagonal dotted lines show the Coulomb 
condition given in Eq. ([3}. 
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FIG. 4: Trajectories of selected non-sliding contacts in the 
system shown in Fig. [T] 



corners of these trapezoids correspond to the times when 
the sliding contact begins or stops sliding. When all con- 
tacts remain non-sliding, the trajectories are no longer 
trapezoids, but straight lines: under loading, each con- 
tact force moves on a straight line, and under unloading, 
it simply retraces its path. The reason for this is given 
in the appendix. 



C. Ratcheting without sliding contacts 

When the force on both walls is varied cyclically: 

= Ly [Po + (t)] , Fy^L:, [Po + Qy (t)] , (4) 

ratcheting can occur without sliding contacts. We car- 
ried out simulations of a small (16 particles) system with 
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FIG. 5: Strain per cycle under elliptic cyclic loading [see 
Eqs. Q and ([Sjl]. The circles are the observed points, and 
the line is a fit of the form A7 — A sin (p. The simulations 
were done with 16 particles. Sliding contacts were surpressed 
by setting /i = 00. 
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FIG. 6: Trajectory of a contact in the {Fn,Ft) plane under 
elliptic cyclic loading, for different values of the phase shift (j) 
(given in degrees). 



Eq. (g]), with 
qx{t) — Aa{l — cos ujt+(i 



qy{t) — A(t(1 — coso^i). (5) 



Note the presence of a phase shift 4> between F^ and Fy. 
We call this form of loading elliptic cyclic loading, be- 
cause ellipses are traced out in the {Fx,Fy) plane. Dur- 
ing these simulations, sliding was suppressed by setting 
fi — 00. The results are shown in Fig. [5l The strain A7 
per cycle is proportional to sin (j). If one traces out the 
path of Qxit), qy{t) in the qx,qy plane, then one obtains 
an ellipse whose area is proportional to sine/). Tracing 
out any contact in the F„ , Ft plane also yields an ellipse 
proportional to sin (j). Some examples are shown in Fig. [S] 
This suggests that ratcheting is related to the area en- 
closed by trajectories in the (F„, Ft) plane. 
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D. Sign of the strain 

Ratcheting with smaU numbers of particles has another 
distinguishing property: the strain accumulation can be 
either positive or negative. Note that in Eq. ([T]), the 
average imposed force does not correspond to an isotropic 
pressure, because q{t) > 0. The pressure exerted by the 
walls on the top and bottom of the sample are larger 
than at the side walls. Thus one expects the sample 
to be gradually flattened, with the top and bottom walls 
moving toward each other, while the side walls are pushed 
apart. This corresponds to 7 < in Eq. A series 
of 100 different ratcheting simulations with 16 particles 
were performed, differing from each other only in the 
initial condition. Of these 100 simulations, 71 exhibited 
unambiguously ratcheting. Of these 71 cases ratcheting, 
51 had 7 < as expected, but the remaining 20 had 
7 > 0. 

On the other hand, when the sample size is larger, one 
has 7 < whenever there is ratcheting. A second series 
of 25 simulations, this time with 400 particles, yields 18 
unambiguously ratcheting simulations, all with 7 < 0. 
The range of 7 observed is also much smaller than for 
16 particles. These results suggest that the strain accu- 
mulated by a large sample is some kind of average over 
strain accumulated by the small regions composing it. In 
these small regions, there can be either negative or pos- 
itive strain, but after averaging, these fluctuations are 
smoothed out, so a large sample has a quite predictable 
behavior. 

Previous studies of granular ratcheting always consid- 
ered large numbers of particles, so this was never noticed. 



III. ORIGIN OF RATCHETING 

A. Particle interaction model 

We now turn our attention from the description of 
granular ratcheting to its cause. We will consider a non- 
sliding contact between two particles subjected to cyclic 
external forces. To facilitate the analysis, we assume that 
Fn is linear in the overlap distance. One imagines that 
when two grains flrst touch, two springs are created, one 
in the tangential and the other in the normal direction. 
Both springs obey Hooke's law so that the normal and 
tangential contact forces are proportional to the spring 
elongations Z)„, Dt: 



Ft 



(6) 



where if„ and Kt are the spring constants. Here, F„ > 
is interpreted as pushing the particles apart, and -D„ < 
occurs when the particles overlap. Eq. ^ holds only for 
touching Dn < particles, so F„ > in accord with 
Eq. ([3]). On the other hand, Dt can have either sign, 
corresponding to the two opposite tangential directions 
(up and down in Fig. [7]) . 




FIG. 7: Definitions for particle interaction laws. Tlie unit 
vector n defined in Eq. (|10|l points from particle j toward 
particle i, and t = z x n, where the z axis points out of the 
page. 



The springs are stretched by the relative motion of 
the particles, as long this does not violate any of the 
conditions in Eq. ([3]). When the contact is non-sliding, 
one has 



dD„ 



dDt 



Vt, 



(7) 



where Vn and Vt are just the relative velocities at the 
point of contact: 



Vn = (v, - Vj) 

Vt = (v.-v,) 



n, 

t - riUJi - rjUj, 



(8) 
(9) 



where v^, tOi, and are the velocity, angular velocity 
and radius of particle i, and i and j label the touching 
particles. The unit vector n 



n 



(10) 



points from particle j toward particle i, and t is a tangent 
vector. If the two-dimensional space is assumed to be 
embedded in a three dimensional one, t can be deflned 
as t = z X n, as shown in Fig. [T] The forces F„ and Ft 
are then directed along n and t respectively. Note that 
the signs in Eq. ^ depend on the choice of n, t, and the 
meaning of positive and negative Dt- In Fig. [71 Dt > 
means points attached to particle i move upward relative 
to points attached to particle j. 

If a contact opens, then Fn ~ Ft ~ in accord with 
the first condition of Eq. ([3]). If two separated particles 
come together again, there is no memory of the previous 
contact. 

The second condition in Eq. Q is enforced by setting 



Dt 



^tJ'-rrDn 



(11) 



5 



whenever using Eq. ([7]) would lead to a violation of 
Eq. ©. 

Sliding contacts are accounted for by modifying 
Eq. jll), but we do not need to consider this in detail, 
since sliding is not needed for ratcheting to occur, as 
shown in Sec. Ill CI 

Note that no damping has been included in Eq. ([6|. 
This is because ratcheting is a quasi-static phenomenon. 
As the frequency of the cyclic loading becomes very long, 
the deformation per cycle approaches a constant. In the 
limit of an infinitely long cycle, the particle velocities 
vanish. Any damping will also vanish, since it is pro- 
portional to the velocities. Since ratcheting exists in the 
limit of infinitely long cycles, one does not need to con- 
sider damping in order to understand granular ratchet- 
ing. Damping is always included in simulations to model 
the loss of energy when grains collide or slide against one 
another. 

The model that has been described above has been in 
use for almost thirty years [13] ■ It has been used in many 
different studies, and considered to be well understood. 
Nevertheless, we show that this model contains an ap- 
proximation that generates granular ratcheting. 



B. Path dependent potential energy 

Granular ratcheting occurs because the model de- 
scribed in Sec. (jlll Ap yields a path-dependent potential 
energy. Here, we are referring to the potential energy 
stored in a contact when two particles overlap. It mod- 
els the elastic energy stored when two grains are pushed 
together. If the force pushing two particles together is 
suddenly released, this potential energy is converted into 
kinetic energy, and the particles will separate. When 
they separate, the highest possible kinetic energy they 
can attain is 

E=^{Kr.Dl + KtD',). (12) 

Thus Eq. (|12p gives the potential energy stored in the 
contact. 

We now show that this energy can be changed if the 
particles execute a closed path relative to one another. 
Consider the path shown in Fig. [51 This figure shows 
a single contact between two particles. Let the lower 
particle be fixed, and let the contact between the two 
grains be always non-sliding. The point A marks the 
center of the upper particle, which is then moved so that 
it traces out the path: A^B^C^D^A. Neither 
particle rotates. Even though the path is closed, the 
length Dt of the tangential spring is changed. 

The changes in Z?„ and Dt during this cycle are 
sketched in Fig. [9l The segments AB and CD change 
only the normal spring length Z)„, whereas the arcs BC 
and DA change only the tangential spring length Dt- 
Segments AB and CD are of equal length, so at the end 




FIG. 8: The origin of granular ratcheting. This figure shows 
two touching disks. The lower disk is fixed, and the upper 
disk moves without rotating, its center tracing out the closed 
path y4— >_B— >C— >-D— >A. The contact forces return to 
their initial state only if the upper particle stops at A' instead 
of proceeding to A. 



of the cycle, Z3„ has returned to its initial value. How- 
ever, arc BC is shorter than arc DA because it lies closer 
to the center of the lower particle. Therefore, Dt does not 
return to its original value, because Eq. ([9]) implies that 
the change in the tangential spring length depends only 
on the distance moved, irrespective of the distance be- 
tween the touching particles. Thus a cycle that returns 
the particles to their initial positions can modify the po- 
tential energy. The potential energy of a contact does not 
depend only on the coordinates of the grains, but also on 
the past relative movements. 

To see why this leads to granular ratcheting, note that 
Dt determines not only the potential energy, but also 
the tangential force. Thus, when the particle executes 
the cycle shown in Fig. [H and returns to A, the contact 
force has also been modified. 

Now let us consider a packing of particles, subjected to 
quasi-static cyclic loading. At the beginning of the load- 
ing cycle, the packing is in static equilibrium, so that the 
net force on each particle vanishes. As the external load 
is varied, the contact forces and the particle positions 
must also change. After one loading cycle, the external 
load has returned to its initial value. If all the particles 
return to their initial positions and all the contact forces 
to their initial values, then there is no deformation of the 
sample, and thus no ratcheting. On the other hand, if 
the contact forces have not returned to their initial val- 
ues, the packing will no longer be in force equilibrium, 
and some deformation must occur. 
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consider a more general motion, such as the one shown 
in Fig. [ini where the upper particle traces out a con- 
vex polygon. The trajectories of the non-sliding contacts 
shown in Fig. [Jare possible examples. Again let us use 
polar coordinates, with the origin placed at O. The path 
of the upper particle will now be given by r{t) and 9{t), 
where t is time. We identify two points, labeled A and B 
in the figure, where 9{t) attains its maximum and mini- 
mum values. At any time, the tangential velocity is given 
by 



Vt 



'It' 



(13) 



and so the total change of the tangential spring, as the 
particle moves from A to i? is 



FIG. 9: Spring lengths D„, Dt for the cycle sketched in Fig. [8] 
Initially, Dt = and Dn < 0. As the upper disk moves from 
A to B, Dn decreases. Then Dt increases as the upper disk 
moves from B to C. When it arrives at D, Fn has returned to 
its original value. However, upon closing the cycle (returning 
to A), Ft does not return to its initial value. 




FIG. 10: A cycle where the lower particle is fixed and the 
upper particle traces out a convex polygon. If the trajectory 
is given in polar coordinates with the origin being the center 
O of the lower particle, points A and B are where the angular 
coordinate takes on its maximum and minimum values. 



C. The role of sliding contacts 

The explanation of ratcheting presented here makes no 
reference to sliding contacts. Yet earlier studies identi- 
fied sliding contacts as a requirement for ratcheting. To 
understand the role of sliding contacts, it is necessary to 



dO /■ ^ 

ADtiA^B)= I r{t)-dt= rAB{e)de, {U) 



where r^_B(0) gives the trajectory that the particle fol- 
lows from j4 to -B. The change in Dt on the return trip 
is 



AA(B 



A)= : 



rBA{e)de, 



(15) 



where rsAid) is the path followed from B back to A. The 
total change in length of the tangential spring is obtained 
by adding Eqs. (fH]) and ([TS]) together: 



[rAB[e)-rBA{e)]de. 



(16) 



Ba 



If the particles are very stiff, then the deformations are 



small: 16*^1 — 9b\ and 
Eq. fTB]) can be written 

AA = 



Vab - rBA\ < r-i -I- r,. Then 



(17) 



where a is the area enclosed by the trajectory of upper 
particle. 

Now the role of the sliding contacts becomes clear. If 
there are no sliding contacts, then the trajectories are 
straight lines, and rAB(O) = rBA{S) for all Ob < < Oa, 
and a = 0. Thus there is no change in Dt if the particles 
return to their original positions, and ratcheting does 
not occur. On the other hand, the presence of sliding 
contacts guarantees that a 7^ 0, so the particles cannot 
return to their original positions, and ratcheting occurs. 
The reason why sliding contacts are required to obtain 
trajectories that enclose a non-zero area is explained in 
Sec. El 



FIG. 11: The definition of the tangential spring. Its length is 
equal to the arc length AC plus the arc length DB. Points A 
and B are defined when the two particles first touch. They 
are carried by the rigid body motion of the particles. Points 
C and D are defined by the intersection of the the line con- 
necting the centers (horizontal line) with the particle surfaces. 



IV. ANGULAR MOLECULAR DYNAMICS 
A. Algorithm 

1. Definition of the tangential spring 

To confirm our explanation of granular ratcheting, we 
show how it can be eliminated by using a new method 
of calculating the tangential forces where the potential 
energy is path- independent. To do so, we retain Eq. (|12p . 
but define Z3„ and Dt in such a way that they depend 
only on the coordinates of the particles. For the spring 
in the normal direction, this is straightforward. If the 
particle positions are given, the overlapping distance can 
be used as the normal spring length: 

Dn = ri + rj - |xi - Xj I . (18) 

where x^ and Xj are the positions of the touching particles 
and Ti and rj their radii. 

For the tangential spring, the point of first contact 
must be stored. Let us imagine that when two particles 
first touch, a spot is painted on each particle, marking 
the point where they touch. Let these points be called 
A and B. The points of first contact are fixed to the 
particle surfaces, and thus carried with the subsequent 
solid-body motion of the particles. To determine the tan- 
gential spring length at a later time, one first determines 
the current points of contact C and D. These points are 
defined by the intersection of the particle surfaces with 
the line connecting the centers. The tangential spring 
length is the length of the arc AC, plus the length of the 
arc DB, as shown in Fig. [TT] 



FIG. 12: Definition of rolling distance. It is the average of the 
distance AC and the distance BD. The points are defined as 
in Fig. [TT] 

One useful side effect of calculating the tangential 
springs in this way is that one can easily obtain the dis- 
tance the particles roll relative to one another. If two 
particles touch, and then roll without sliding, the points 
A through D will be as sketched in Fig.[T21 The distance 
rolled is the length of the arc AC or BD (see Fig. [T2|) . 
Note that this measure of the rolling is objective, because 
it is based on points fixed on the particles themselves, and 
is thus independent of any solid-body motion imposed on 
the two particles. 

Each point of a circle can be assigned an angle, ob- 
tained by measuring the angle between the cc-axis and a 
line determined by the point in question and the center 
of the circle. In this way one can assign the angles Oa, 
9b, Oc, and 9d to the points A, B, C, D, and arc lengths 
can now be calculated by subtracting angles. Thus the 
arc length AC is r.i{6A — Oc)- Note that the arc length 
has a sign, which is necessary for distinguishing between 
rolling and sliding. Note also that DB is rj{6B — Od) 
since angles are measured with respect to particle j and 
not i. 

Now we can write the tangential spring length as 

Dt = r,{9c - Oa) + n{0D - Bb), (19) 

where we adopt the convention that Dt increasing means 
that point A in Fig. llll moves upward (i.e., 9a decreases), 
and B moves downward (i.e., Ob decreases). The rolling 
distance is 

„ r,{0c-9A)-n{9D-0B) 

-i^roU = ^ \M) 

2. Direct implementation 

The most obvious way to implement this algorithm 
is to calculate the angles 9a through 9d, and then use 



8 



Eq. (fT9)) to calculate the spring length. 

The angles 9a and 6b can be found by integrating the 
equations 9a — tOj, 6b = oji. But it may be more econom- 
ical to assign an angular coordinate (t) to each particle 
i. When the particle positions are updated, 9i can be 
updated as well, using 9i — oji. Then at the time of 
first contact, one stores 6i{t^), and A9i = 6i{t^) — 6'„, 
where 0„ is the angle of the point of first contact at time 
t*. Then at any later time t: 



6A = 6,{t)-9,{U) + A6, 



(21) 



The angles 6c' and 6d are calculated at each time step 
from the positions of the particles. Writing and Uy 
for the two components of n. 



9c ^ 



tan 1 n^/ny, > 0, 

TT — tan~^ Ux/riy Ux < 0. 



(22) 



Then 6c is moved into the correct interval by adding or 
subtracting 2tt. Then one uses 9]j = 9c ^ tt. In this way, 
both the rotation and translation of the particles is taken 
into account. 

If a contact slides, one moves the points A and B along 
the particle surfaces so that Eq. pTjl is satisfied. In a 
similar way, one could move these points to set DroU = 
while leaving Dt unchanged. 



3. Implementation through integration 

An alternative way of implementing this algorithm is 
to modify the existing Cundall-Strack algorithm. A few 
minor modifications are necessary to obtain an algorithm 
with a potential energy given in Eq. (|12p which is path- 
independent, with Dt defined as in Eq. 

To do this, first express Dt and Dro\i in terms of the 
motion of the particles. To do this, we need 



n 



(V, - Vj) -t 



(23) 



where the tangent vector is t = z x n as in Fig. [T] Then 
Eqs. (dH), (EI]), and (HH) give 

Dt = -riUJi - TjUjj + a{-Vi - Vj) ■ t, (24) 

where we have defined 

a = ^1±^ (25) 



Note that the first of these is equivalent to Eq. ([7|) and ([9|) 
only when a = 1 or jx^ — x-, | = + rj^ i.e. when the par- 
ticles are just touching. The usual implementation of the 
Cundall and Strack model, therefore, contains an approx- 
imation, namely a « 1. Normally one chooses a stiffness 
high enough so that this approximation is reasonable, but 
it nevertheless has an effect on the simulation results. 



In the same way, one can obtain a rolling velocity from 
Eq. (1201): 



-Droll 



(26) 



To obtain the equations of motion, one cannot simply 
use Eq. To guarantee conservation of energy, one 

defines the Lagrangian Q 



L = T^V 



(27) 



where T is the kinetic energy of a system, and V is the 
potential energy. In our case, we consider the two touch- 
ing particles whose kinetic energy is 



1 



1 



1 



1 



T = -TO, (±2 + yf) + -TO, (.t2 + 2)2) + _I^02 + _I^02 

(28) 

The potential energy V is given by Eq. (fT^ . The equa- 
tions of motion are then given by 



A - — - 

dt \dq J dq 



(29) 



where q is one of the coordinates of the grains 
Xi,yi,9i, Xj,yj, 9j. Applying this equation yields 



= ±KnDnn ± aKtDtt, 



= KtDtr, 



(30) 



Note that these differ from Eq. ([6]) by the presence of 
the factor a in the tangential force. This same factor 
appears in Eq. (p4|) . This means that this new method 
can be implemented simply by inserting this factor in the 
appropriate places in the program. 



B. Results 

We have compared the traditional Cundall-Strack al- 
gorithm used in Sees. Ill CI and IIIDI with the the two 
different implementations of the angle-based algorithm 
discussed in Sec. IIV Al 



1. Simulation Parameters 

In all cases, the initial condition was generated by plac- 
ing grains on a lattice in a square domain. The radii 
are uniformly distributed within the interval rniax[0.7, 1], 
where Tmax is chosen so that the desired number of par- 
ticles will fit in the domain. 

Two walls of the domain are fixed, and the other two 
are movable. A force proportional to wall length is ap- 
plied to the movable walls, and they compress the par- 
ticles at uniform stress into a packing. During this time 
of compression, the particles are smooth (friction ratio 
fi = 0). Once this compression is complete, one sets 
/i = 0.2, and imposes cyclic loading as described above. 
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FIG. 13: Deformation per cycle for all ratcheting simula- 
tions. 100 different configurations were generated using dif- 
ferent random number seeds. 

The system of units for the simulation is given by the 
initial length L of the system, the (two-dimensional) pres- 
sure p applied during compression, and the density p of 
the particles. In these units, the stiffness of th e par ticles 
is Kn — Kt — lOOp. The unit of time is r — ^Jm/p. One 
cycle lasts lOr. At least 100 cycles were performed in all 
simulations. 

The position of the movable walls are recorded at the 
time of minimum force during each cycle. By comparing 
these values from one cycle to the next, an accumulation 
of strain can be detected. To determine whether a sam- 
ple ratchets, the following procedure was applied. First, 
the first 29 cycles were neglected to eliminate transients. 
Then L^o and LyQ were defined by the positions of the 
walls at the beginning of the thirtieth cycle. Next, the 
strain 7, defined in Eq. ([2]) was calculated for each sub- 
sequent cycle. Finally, we checked whether 7 increases 
linearly with cycle number N . This was done by fitting 
a line to the observed (7,iV), and calculating the root 
mean square deviation of the observed points from the 
fit. If this number was smaller than the slope, the sim- 
ulation was judged to exhibit ratcheting. Otherwise, it 
was considered non-ratcheting. 

2. Ratcheting in small systems 

We subjected 100 different small packings (16 parti- 
cles) to cyclic loading as described above. With the 
unmodified Cundall-Strack algorithm, 71 simulations ex- 
hibited ratcheting. The deformation per cycle A7 varied 
over a wide range: 10^^^ < IA7I < 10~^, with a geomet- 
ric mean of 8 x 10^^. Both positive and negative values 
were observed: — 10~^ < A7 < 4 x lO""". 

When the Cundall-Strack algorithm is modified as de- 
scribed in Sec. IIV A3[ 62 still simulations exhibit ratch- 
eting, but at a much lower amplitude. One observes 
10~^^ < IA7I < 2.5 X 10""'^° with a geometric mean of 
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FIG. 14: Dependence of ratcheting on time step. Crosses: 
original algorithm. Circles: corrected algorithm. The time 
step is given in multiples of ^ K/m-mi-n, where mmin is the 
smallest particle in the simulation. 

1.1 X 10^^^. These results are summarized in Fig. [T31 
One sees that the use of the corrected equations leads to 
a 10'*-fold reduction of granular ratcheting. 

The remaining granular ratcheting is due to integra- 
tion errors. This can be shown by taking a single initial 
condition, and changing the time step. Typical results 
are shown in Fig. [TH With the original Cundall-Strack 
algorithm, ratcheting is independent of the time step. 
When it is modified, then the ratcheting deformation is 
proportional to the time step. 

Finally, when the method described in Sec. IIV Al ls im- 
plemented by direct calculation of angles, no simulations 
ratchet. 28 of the simulations exhibit a constant strain 
with I7I < 10"^** for every cycle. Note that such small de- 
formations are not even visible on Figs. [T51 and [HI The 
others exhibit a variety of other behaviors that will be 
discussed in the next section. 



3. New behaviors in small systems 

Once ratcheting has been removed or reduced, new 
behaviors come to the foreground. The most common 
of these are outliers. The strain is independent of cy- 
cle number, except for occasional cycles. An example is 
shown in Fig. [151 

A closer inspection of these simulations reveals that 
these outliers are due to "rattlers" : particles without con- 
tacts. Since there is no gravity, rattlers float inside cages 
in the packing. Occasionally, they collide with walls of 
their cage. These collisions coincide with the outlying 
points. Once the collision is past, the packing returns to 
its initial state, and the rattler floats off toward another 
part of its cage. The packing that produced Fig. [T5| is 
shown in Fig. 1161 The rattler is found in the upper cen- 
ter of the packing. It moves within a cage formed by 
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FIG. 15: Strain as a function of cycle number for a simulation FIG. 17: Accumulated strain 7 versus cycle number for a sim- 
with outliers. ulation performed with the angle based method of calculating 

tangential forces. The strain is periodic. This arises when the 
rattler is confined in a very small cage. 



In other cases, every point can be considered as an out- 
lier. Another thing that can happen when the rattler's 
cage is small is that it can be driven in the cage by the 
motion of the surrouding particles in a periodic way. This 
leads to a periodic dependence of 7 on iV, as shown in 
Fig. [171 

Rattler-induced outliers exist also in the original Cun- 
dall and Strack method. When on inspects the 29 non- 
ratcheting simulations, one finds that 25 of them have 
perturbations due to rattlers. 

Another effect that rattlers can cause is a sudden step 
in the strain. This occurs when the particles forming 
the cage have only weak contact forces. The rattler can 
induce a sudden step in the strain by provoking a small 
re-arrangement of these particles. 

One way to reduce the effect of the rattlers is to apply a 
weak gravitational field so that they can no longer float 
slowly from one side of their cage to another. When 
this is done, outliers still exist, but the perturbations 
they introduce are much smaller - 0(10^^^) instead of 
0(10-7). 

Another phenomenon is "shakedown" , which has al- 
ready been investigated [3| • In "shakedown" , the accu- 
mulated strain per cycle decreases each cycle. In FiglT^ 
we show an example. The time required to reach a level 
of negligible strain accumulation is variable. In Fig. 1181 
strain is still accumulating, even after 1000 cycles. In 
other simulations, shakedown occurs after very few sim- 
ulations. 




FIG. 16: A configuration that causes outliers. The arrow 
shows the motion of the rattler. The normal forces are shown 
by the lines connecting particle centers. The tangential forces 
are shown by lines perpendicular to the normal forces. Note 
that the rattler does not participate in the force network. 



five particles and the upper wall. The outlying points in 
Fig.HHcoincide with collisions between the rattler and its 
cage. However, not all collisions leave a trace in Fig. [TSl 
the amplitude registered in Fig. [15] probably strongly de- 
pends on the time within the cycle where the collision 
takes place. 

The frequency of collisions with the cage can vary 
widely. Sometimes only one point out of 70 is perturbed. 



4- Large systems 

To study these phenomena in larger systems, a se- 
ries of 25 simulations with 400 particles was carried out. 
With the original Cundall-Strack method, ratcheting is 
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FIG. 18: Accumulated strain 7 versus cycle number for a sim- 
ulation performed with the angle based method of calculating 
tangential forces. This sample exhibits shake down. 

observed, with a much narrower range of strain accumu- 
lation, as discussed in Sec. Ill Dl When either of the mod- 
ifications proposed in Sec. II V Al ls used, outliers dominate 
the stress - cycle graphs. This is probably because as the 
packing becomes large, the probability of having a rat- 
tler approaches unity. When a weak gravitational field is 
applied, the phenomena of "shakedown" dominates. 

V. CONCLUSIONS 

We have uncovered the cause of granular ratcheting. 
It is due to a potential energy that depends not only on 
the particle positions, but also on their past trajectories. 
As a granular assembly is subjected to cyclic loading, it 
is impossible to return both the particle positions and 
the contact forces to their initial values, so a small defor- 
mation occurs with each cycle. It follows that granular 
ratcheting can be eliminated by defining a potential en- 
ergy depending only on the current particle positions. 
This was confirmed by extensive numerical simulations 
using two different implementations of this idea. 

This result suggests that contact modeling should fo- 
cus on the potential energy as the fundamental quantity, 
and then use Eq. ([M)) to obtain the forces. In contrast, 
the most common approach taken in the literature is to 
directly postulate forces based on physical grounds, with- 
out considering the potential energy. 

One possible criticism of this work is that it is only 
concerned with disks, whereas ratcheting has been found 



in packings of polygons. The motion of disks is much 
simper to analyze, because the rotation of the disks is 
uncoupled from the translation. In polygons, this is no 
longer true. But all that is required is that the poten- 
tial energy be path dependent. One common force law, 
used in granular ratcheting studies @ is to assume that 
the force between polygons is proportional to the overlap 
area. This force law has been shown to violate energy 
conservation Thus it seems likely that the explana- 
tion of granular ratcheting presented here also applies to 
ratcheting of packings of polygons. 

At a detailed level, the numerical mechanism cannot be 
the same as the physical one. Numerical granular ratch- 
eting is a consequence of the way the tangential spring 
is stretched. In experiments, the contacts between the 
particles are not governed by the stretching of springs. 
Indeed, if two touching particles can be considered as 
making up a single elastic body, then force and position 
cycles will coincide, as there is a potential energy. 

However, the results of this paper do show that gran- 
ular ratcheting in the experiments will occur if force and 
position cycles are not equal. In principle, this could be 
checked by examining a single contact under cyclic load- 
ing. Such an experiment would be difficult to do, since 
very small relative displacements must be measured. And 
it must also be mentioned that the idea of two contact- 
ing particles acting as if they were welded together at 
the contact surface is itself an idealization. There may 
be zones of slip at the contact (even when the contact as 
a whole does not slide), and this may give rise to a com- 
plicated behavior when the contact is subjected to cyclic 
loading. Another possibility is that fluid could coat the 
surfaces of the touching particles, possibly lubricating 
them. Or abrasion at the contact point could generate 
very tiny particles trapped between the two touching sur- 
faces. These particles could act like fault gouge [7] on a 
very small scale, facilitating a relative tangential motion. 
All of these effects may lead to a history dependent po- 
tential energy, and thus to granular ratcheting through 
the mechanism discussed in this paper. 
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APPENDIX A: STIFFNESS MATRIX THEORY 

The section presents a very brief review of stiffness 
matrix theory. This theory applies to granular pack- 
ings under quasi-static loading, and thus is applicable to 
granular ratcheting. We explain below how this theory 
explains certain key properties of packing under cyclic 
loading, namely, 

• why particle trajectories are straight lines when 
there are no sliding contacts, and the forcing is 
given by Eq. IJ), 

• why this is no longer true when there are sliding 
contacts, and 

• why the forcing given in Eq. (|4]) generates particle 
trajectories with a non-zero area. 



1. Introduction to stiffness matrix theory 

In stiffness matrix theory , the behavior of the pack- 
ing is piece-wise linear. Thus time can be divided into 
intervals [ti,ti^i] during which the velocities of the par- 
ticles are linearly related to the change in forces: 



dt 



= kv, 



(Al) 



where f^xt represents the external forces {Fr^ and Fy for 
the biaxial box), v contains the velocities of the particles 
and walls, and k is called the stiffness matrix. It relates 
the velocities (or displacement increments) of the parti- 
cles to the change in the force exerted on each particle 
by its neighbors. 

The motion is only piece-wise linear because the stiff- 
ness matrix k depends on the status (sliding or not) of 
each contact. Whenever a contact status changes, k also 
changes. Therefore, the times {ti} which define the inter- 
vals of linearity are the times when one or more contacts 
change status. 

Eq. (jAip holds when the forcing is quasi-static, and 
the particles are quasi-rigid. Quasi-static forcing means 



that the time scale associated with the forcing is much 
longer than the time the packing needs to react. Particles 
can be said to be quasi-rigid if their stiffness is much 
greater than the confining pressure. Note that these two 
assumptions are related: if the particles are very stiff, the 
speed of sound is very high, and the packing can quickly 
react to changes in the external load. 



2. Application to biaxial test 

If one considers a biaxial test, with the forcing given 
by Eq. ([T]), then only the entries of foxt corresponding 
to the walls are non-zero, because no external forces are 
applied to the particles. Furthermore, in Eq. (jAl[) . only 
those components associated with varying forces survive 
differentiation by time. Thus Eq. (|A1|) becomes 



dq 



kv, 



(A2) 



where lix is a vector, all of whose components are zero, 
except the y component of the force on the upper and 
lower walls. All other components of foxt are either zero 
or constant. 

One would like to invert k and bring it onto the left 
hand side of the equation. But k is singular because 
there are certain collective motions that do not change 
the spring lengths, and hence do not change the forces. 
One example is the uniform motion of all particles. They 
do not move relative to one other, and provoke no change 
in force. Let B be the set of all such motions. We can be 
sure that the left hand side of Eq. (jA2[) is orthogonal to 
every member of B, for if it were not, the packing would 
be unstable [1]. 

Now define the matrix k that will act like an inverse 
of k. It is defined by 



kx = f and x _L 



kf 



(A3) 



This equation gives the result of applying k for 3A^ — 
dimB linearly independent vectors. To fully determine 
k, we must say how it acts on the other dim B dimensions 
of R3^. Let F be the range of k. Then: 



kf = 0, for { ±¥. 



(A4) 



This determines k. Note that kk is a projector that 
removes B. 

Using this in Eq. (|A2|) . one can write: 



dq- 

which can be integrated: 

X = xo + g(t)kL:i 



(A5) 



(A6) 



Here x is a vector containing the positions of the parti- 
cles and walls. Eq. (|A6I shows that the position of each 
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particle moves back and forth on a line defined by the 
appropriate components of kL^. . So no area is traced out 
by position cycles, and ratcheting does not occur. 

But if there are sliding contacts, then k does not re- 
main constant. Recall that k changes whenever a contact 
changes status. Thus when a contact starts or stops slid- 
ing, the relation between v and foxt changes, and the 
particle motion changes direction. This is what we saw 
in Fig. m 

Another way to get paths that are not lines is add 
another term to the left hand side of Eq. (|A2p . When we 
use the forcing given in Eq. ^ and (O, then Eq. (|A2[) 
becomes 



and thus the motion is 



^L. + = kv. 



dt 



(A7) 



X = xo -I- ga;(t)kLa; -I- qy{t)'kLy 



(A8) 



Thus if the forcing traces out an area in the {qx,qy) plane, 
than each contact traces out a proportional area in the 
relative position plane. This, together with Eq. [171 ex- 
plains the result in Fig. [5l 



